/**
 * Trabalho 3 de Programacao Concorrente.
 * Prof. Dr. Julio Estrella
 * 
 * Resolucao de sistemas lineares.
 * 
 * Grupo 07 - Turma B
 * Fabio de Albuquerque Dela Antonio 6792558
 * Leonardo Henrique Omoto 6793160
 * Pedro Paulo Canto Martucci 6793100
 * 
 * main.c
 */

#include <geral.h>
#include <matrix.h>
#include <jacobi.h>

#ifdef SEQUENCIAL

/* Versao Sequencial */

int main(int argc, char *argv[]) {

    /* Confere entrada. */
    if(argc < 2) {
        fprintf(stderr, "Uso: %s <arquivo do sistema>\n", argv[0]);
        return 1;
    }
    
    /* Carrega sistema na memoria */
    Sist * sistema = carregaSistema(argv[1]);
    
    int i;
    int n_iteracoes;
    
    double * k_anterior = malloc(sistema->J_ORDER*sizeof(double));
    double * k;
    double erro;
    double aux;
    
    if(k_anterior == NULL) {
        fprintf(stderr, "Impossivel alocar vetor de resultados.\n");
        return 1;
    }
    
    /* Inicia solucoes iniciais com zero. */
    for(i = 0; i < sistema->J_ORDER; i++)
        k_anterior[i] = 0.0f;
    
    /* Faz iteracoes do Jacobi. */
    for(i = 0; i < sistema->J_ITE_MAX; i++) {
    
        k = (double *)malloc(sistema->J_ORDER*sizeof(double));
    
        if(k == NULL) {
            fprintf(stderr, "Impossivel alocar memoria para variaveis.\n");
            return 1;
        }
    
        jacobiRichardson(sistema, k, k_anterior, 0, sistema->J_ORDER);
        
        if((erro = erroJacobi(k, k_anterior, sistema->J_ORDER)) < sistema->J_ERROR) {
            free(k_anterior);
            break;
        } else {
            free(k_anterior);
            k_anterior = k;
        }
    }
    
    /* Armazena numero de iteracoes. */
    n_iteracoes = i;
    
    /* Confere solucao com a posicao dada (J_ROW_TEST) */
    aux = 0.0f;
    
    for(i = 0; i < sistema->J_ORDER; i++)
        aux += sistema->matrix[sistema->J_ROW_TEST][i] * k[i];
    
    /* Imprime a saida (apenas se o valor estiver correto!) */
    //if((aux <= sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER] + erro) && (aux >= sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER] - erro)) {
        printf("-------------------------------------------------------\n");
        printf("Iterations: %d\n", n_iteracoes);
        printf("RowTest: %d => [%.6lf] =? %.6lf\n", sistema->J_ROW_TEST, aux, sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER]);
        printf("-------------------------------------------------------\n");
    //}
    
    /* Desaloca sistema da memoria */
    limpaSistema(sistema);
    
    /* Desaloca memoria do vetor de solucoes */
    free(k);
    
    return 0;
}

#else

#include <mpi.h>
#include <omp.h>

/* Versao Paralela com OpenMP e MPI */

int main(int argc, char *argv[]) {

    int rank, size;
    
    /* Inicializacao da MPI */
    int mpicheck = MPI_Init(&argc, &argv);
    
    if(mpicheck != MPI_SUCCESS)
        return 1;
    
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    
    /* Confere entrada. */
    if(argc < 2) {
        fprintf(stderr, "Uso: %s <arquivo do sistema>\n", argv[0]);
        MPI_Finalize();
        return 1;
    }
    
    /* Carrega sistema na memoria */
    Sist * sistema = carregaSistema(argv[1]);
    
    /* Variaveis do problema */
    double * k_anterior = malloc(sistema->J_ORDER*sizeof(double));
    double * k;
    int i;
    
    /* Variaveis do paralelismo */
    MPI_Status status;
    int tag = 999;
    int continua = 1;
    int inicio, fim;
    
    /* Verifica alocacao de vetor de solucoes anterior. */
    if(k_anterior == NULL) {
        fprintf(stderr, "Impossivel alocar vetor de resultados.\n");
        MPI_Finalize();
        return 1;        
    }
    
    /* Execucao do no mestre. */
    if(rank == MASTER) {

        #ifdef DEBUG
        printf("Mestre!\n");
        #endif
        
        int n_iteracoes;
        int j, w;
        
        double erro;
        double aux;
        
        double * k_escravo = malloc(sistema->J_ORDER*sizeof(double));
        k = malloc(sistema->J_ORDER*sizeof(double));
        
        if(k == NULL || k_escravo == NULL) {
            fprintf(stderr, "Impossivel alocar vetor de resultados.\n");
            MPI_Finalize();
            return 1;            
        }
        
        /* Inicia solucoes iniciais com zero. */
        for(i = 0; i < sistema->J_ORDER; i++)
            k[i] = 0.0f;
        
        /* Faz iteracoes do Jacobi. */
        for(i = 0;;i++) {
            
            /* Envia solucoes da iteracao anterior para nos escravos... */
            MPI_Bcast(k, sistema->J_ORDER, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
            
            /* Armazena solucao anterior (para conferir posteriormente). */
            for(w = 0; w < sistema->J_ORDER; w++)
                k_anterior[w] = k[w];
            
            /* Recebe solucoes novas dos escravos e agrega na solucao do mestre. */
            for(j = 1; j < size; j++) {
                
                inicio = (j - 1)*(sistema->J_ORDER/(size - 1));
                fim = (j == size - 1) ? sistema->J_ORDER : (j)*(sistema->J_ORDER/(size - 1));
                
                MPI_Recv(k_escravo, (fim - inicio), MPI_DOUBLE, j, tag, MPI_COMM_WORLD, &status);
                
                for(w = inicio; w < fim; w++)
                    k[w] = k_escravo[w - inicio];    
            }
            
            /* Calcula erro e criterios de parada... */
            if((erro = erroJacobi(k, k_anterior, sistema->J_ORDER)) < sistema->J_ERROR) {
                /* Envia sinal para os escravos terminarem a execucao. */
                continua = 0;
                MPI_Bcast(&continua, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
                
                break;
            } else if(i == (sistema->J_ITE_MAX - 1)) { 
                /* Envia sinal para os escravos terminarem a execucao. */
                continua = 0;
                MPI_Bcast(&continua, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
                
                i++;
                break;
            } else {
                /* Envia sinal para os escravos continuarem a execucao. */
                continua = 1;
                MPI_Bcast(&continua, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
            }
        }
        
        /* Armazena numero de iteracoes. */
        n_iteracoes = i;
        
        /* Confere solucao com a posicao dada (J_ROW_TEST) */
        aux = 0.0f;
        
        for(i = 0; i < sistema->J_ORDER; i++)
            aux += sistema->matrix[sistema->J_ROW_TEST][i] * k[i];
        
        /* Imprime a saida (apenas se o valor estiver correto!) */
        //if((aux <= sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER] + erro) && (aux >= sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER] - erro)) {
            printf("-------------------------------------------------------\n");
            printf("Iterations: %d\n", n_iteracoes);
            printf("RowTest: %d => [%.6lf] =? %.6lf\n", sistema->J_ROW_TEST, aux, sistema->matrix[sistema->J_ROW_TEST][sistema->J_ORDER]);
            printf("-------------------------------------------------------\n");
        //}

        /* Desaloca vetor de solucoes. */
        free(k);
        free(k_escravo);
        
    }
    
    /* Execucao nos demais nos. */
    else {
        
        /* Intervalos de processamento do no escravo. */
        inicio = (rank - 1)*(sistema->J_ORDER/(size - 1));
        fim = (rank == size - 1) ? sistema->J_ORDER : (rank)*(sistema->J_ORDER/(size - 1));
        
        /* Aloca vetor local de solucoes. */
        k = (double *)malloc((fim - inicio)*sizeof(double));
        
        if(k == NULL) {
            fprintf(stderr, "Impossivel alocar vetor local de solucoes.\n");
            MPI_Finalize();
            return 1;
        }
        
        /* Ponteiro auxiliar */
        double * k_aux = &k[(fim + inicio)/2 - inicio];
        
        #ifdef DEBUG
        printf("Escravo! (%d a %d)\n", inicio, fim);
        #endif
        
        /* Loop do no escravo. */
        while(continua) {
            
            /* Recebe solucao anterior do no mestre. */
            MPI_Bcast(k_anterior, sistema->J_ORDER, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
            
            /* Efetua calculo de sua parte. */
            #pragma omp parallel
            {
                /* Utiliza decomposicao de dados local em duas threads. */
                #pragma omp sections 
                {
                    #pragma omp section
                    {
                        jacobiRichardson(sistema, k, k_anterior, inicio, (fim + inicio)/2);
                    }
                    
                    #pragma omp section
                    {
                        jacobiRichardson(sistema, k_aux, k_anterior, (fim + inicio)/2, fim);
                    }
                }
            }
            
            /* Envia sua parte para o no mestre. */
            MPI_Send(k, (fim - inicio), MPI_DOUBLE, MASTER, tag, MPI_COMM_WORLD);
            
            /* Recebe sinal de continuidade do no mestre. */
            MPI_Bcast(&continua, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
            
        }
     
        /* Desaloca vetor de solucoes local */
        free(k);
        
    }

    /* Desaloca vetor auxiliar de solucoes anteriores. */
    free(k_anterior);

    /* Desaloca sistema da memoria */
    limpaSistema(sistema);

    /* Termina MPI */
    MPI_Finalize();
    
    return 0;
}

#endif
